Quantum order by disorder in a spin-one frustrated magnet on the kagome lattice 
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We study the XXZ spin-one quantum magnet on the kagome lattice as an example where quantum 
fluctuations on highly degenerate classical ground states lead to various exotic quantum ground 
states. Previous studies have predicted several quantum phases, but different analytical approaches 
do not necessarily lead to the same physical picture. In this work, we use Quantum Monte Carlo 
computations to critically examine some of the predictions made in the string-net mean-field theory 
and the degenerate perturbation theory combined with duality analysis and effective field theory. 
It is found that the resulting phase diagram differs from some of the previous predictions. Further 
implications of our results to different analytical approaches are discussed. 
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I. INTRODUCTION 

Macroscopic degeneracy of classical ground states in 
frustrated magnets is a fertile ground for emergence of 
unusual quantum ground states such as spiral magnetic 
order, quantum spin liquid, and valence bond solid (VBS) 
that occur via quantum fluctuationSfi Any perturbation 
on the degenerate classical ground states, however, is in- 
herently strong and there may also be more than one 
competing ground states that are extremely close in en- 
ergy. The situation is very reminiscent of quantum Hall 
states that emerge from highly degenerate Landau lev- 
els. This is indeed one of the reasons why identification 
of true quantum ground state in frustrated magnets is 
such a difficult task. 

Different kinds of non-perturbative analytical ap- 
proaches have been proposed in literature to understand 
this quantum order (or disorder) by disorder phenomena. 
These approaches include various degenerate perturba- 
tion theories duality analysis combined with effective 
field theories fii^ and string net condensation picture.^ 
Since different approaches may not necessarily lead to the 
same conclusion, it is important to understand the limi- 
tations of various approaches. One useful way to obtain 
such information is doing unbiased numerics by appropri- 
ately choosing concrete examples where predictions from 
different approaches may vary. 

Here we consider the following spin-one XXZ model on 
the kagome lattice as such an example. 



z\2 



H = -Jj_ Y^iS^SJ + SfS]) + J, 

(1) 

where J±,Jz^D > 0, the first two sums run over the 
nearest neighbors on the kagome lattice, and D is the 
strength of the single-ion anisotropy. This model is also 
equivalent to a boson model with nearest neighbor re- 
pulsive interaction where boson occupation number can 
only assume 0,1,2. Previously two different analytical 
approaches have been used to study the phase diagram 
of this model. Using degenerate perturbation theory and 
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FIG. 1: (color online). The schematic phase diagram. First 
order phase transitions are denoted by thin lines (blue and 
green online) and a continuous phase transition is denoted by 
a thick line (red online). 



mapping to a dimer model combined with duality analy- 
sis, Xu and Moored predict that there exist three differ- 
ent phases: an XY ferromagnetic phase for J± ^ Jz,D, 
a plaquette valence bond solid (VBS-H) phase with res- 
onating hexagons for J± <^ D < Jz, and a gapped "pho- 
ton" (quantum paramagnetic) phase for J± <^ Jz < D. 
The gapped "photon" phase is a descendant of an un- 
stable two-dimensional quantum spin liquid phase with 
linearly dispersing neutral "photon" modes. It was also 
predicted that there would be a direct transition between 
VBS-H and a gapped "photon" phase. Using a quite 
different approach. Levin and Wen® proposed a mean- 
field theory where a class of variational wavefunctions 
based on the so-called "string-net" picture was used to 
map out the global phase diagram. While they also 
predict the existence of a plaquette phase as well as a 
gapped "photon" phase, their plaquette phase is char- 
acterized by frozen spin configurations and this frozen 
plaquette (FP) phase is different from VBS-H in lattice 
and spin symmetriesi^ In this string-net mean-field the- 
ory, the "photon" phase corresponds to the string-net 
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condensed phase^ and the unstable two-dimensional spin 
liquid phase mentioned above. This "photon" phase was 
found in a narrow region around Jz ^ D if one ignored 
the nonperturbative instanton contribution. The "pho- 
ton" excitation acquires a finite gap, however, due to the 
instanton effect in 2+1 dimensions.— The resulting state 
is an ordinary quantum paramagnetic state. 

In this work, we map out the phase diagram of 
the model given by Eq. [T] using quantum Monte Carlo 
method. We use a plaquette generalization^ of the 
stochastic series expansion algorithm.^ We consider only 
the parameter region J± > 0, > 0, and D > 0. The 
schematic phase diagram is shown in Fig. [T] First, we 
find that a plaquette ordered VBS phase indeed arises 
and it is the VBS-H phase, not the FP. In addition to 
the phases discussed in the previous works, we also find 
an additional phase with resonating lattice units for small 
values of D/J^ and J±/ Jz- The ground state in this re- 
gion is either a phase with resonating triangles (VBS-T) 
or a phase with resonating bow ties (VBS-B). It was ar- 
gued in Ref. [sl that there should be a direct continuous 
transition from the VBS-H phase to the paramagnetic 
(gapped "photon") phase. We do not observe such a 
transition for J±/Jz ^ 0.08. However, we cannot fully 
rule out the possibility of a direct VBS-H-paramagnet 
phase transition at even smaller values of J±/ Jz- 

The rest of the paper is organized as follows. In Sec- 
tion II, we will discuss the details of the phase diagram 
and properties of all the phases discovered in the nu- 
merics. Here we also discuss the difference between our 
results and the predictions of previous works. We dis- 
cuss the implications of our results to various analytical 
approaches in Section HI. 



II. CONSTRUCTION OF THE PHASE 
DIAGRAM 




FIG. 2: (color online). The correlation function Cb(ri — rs) 
between the bond indicated by a large text (red online) and 
the other bonds for L = 24, J^/J± = 10, D/J^ = 0.4, and 
T — Jx/48. Resonating hexagons are denoted by thick lines. 
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FIG. 3: (color online). Finite size scaling of the equal time 
bond-bond structure factor S'b(Qo) for L = 24, Jz/J± = 10, 
D/Jz = 0.4, and T = J±/48. In this and the other structure 
factor plots, error bars (if not visible) are smaller than the 
symbol sizes, and the line shows a linear extrapolation to the 
thermodynamic limit. 



A. VBS-H phase 

We begin with the analysis of the VBS-H phase that 
is found in a wide lobe as shown in Fig. [TJ This phase 
has plaquette order. Two different plaquette phases were 
proposed previously for this region of parameters. In the 
plaquette phase (VBS-H) predicted by Xu and Moore^^ 
spins resonate on a third of the total number of hexagons. 
These spins have {S^) = 6/2 around resonating hexagons, 
where b is some number. The rest of the spins is fixed 
to (Sf) = —b. In the frozen plaquette phase (FP) pre- 
dicted by Levin and Wen,^ spins are frozen on a third 
of all the hexagons with alternating values of (S^) = ±a 
around those frozen hexagons; the rest of the spins has 
(5f ) = 0.® Thus, these two phases break Ising and lattice 
symmetries in different ways. Note that these two phases 
can also be distinguished by the distribution function of 
(S'f ) (more discussions are given below). 

In the resonating plaquette phase, the bond-bond cor- 
relation functions should have peaks in momentum space 



corresponding to the specific pattern in real space char- 
acterizing plaquette bond order. In Fig. ^ we show 
the bond-bond correlations function in real space that 
is given by 



i / B^rdr - Bo 



^ / Bsrdr - Bo 



(2) 

where — J±{SfSj + SfSj) is the off-diagonal 

bond operator (at imaginary time r) of the bond a con- 
necting spins i and j, and Bo denotes the background 
bond strength. In Fig. [21 one can clearly see the pat- 
tern that is compatible with the VBS-H phase. To con- 
firm that the bond order is long-ranged, we show in 
Fig. [3] the finite size scaling of the equal time bond- 
bond structure factor S'b(q) at the ordering wavevector 
q=Qo = (47r/3,0), 

^b(q) - N{Bl,B^r), (3) 
where Bqr = {l/N)J2a '2xp{iqra). Note that the 
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FIG. 4: (color online). Finite size scaling of the equal time 
spin-spin structure factor S(Qo) for L — 24, Jz/Ji_ = 10, 
D/Jz = 0.4, and T = Jj_/48. 



FIG. 6: (color online). Finite size scaling of the equal time 
bond-bond structure factor S'b(Q3) for Jz/Jx — 10.5, D = 0, 
and T = Jx/12. 
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FIG. 5: (color online). Probability distribution (arbitrary 
units) of Si for L = 48, Jz/Jx = 10, D/J^ = 0.4, and T = 
Jx/48. The distribution function is independent of the system 
size (not shown). 




expression of B^t involves the sum over the bond index 
a and N is the number of sites. The structure factor di- 
vided by the number of sites clearly goes to a finite value 
in the thermodynamic limit even though the extrapolated 
value is quite small. 

On the other hand, the bond-bond correlation de- 
scribed in the previous paragraph might also be com- 
patible with the FP phase. More robust diagnostic is 
the distribution function of local (S'f), which should 
have two peaks for the VBS-H phase and three peaks 
for the FP phase if the system is in one of the six 
degenerate symmetry-broken states.^- For each Monte 
Carlo configuration, we compute the time-averaged S'f — 

(1//3) drS^^, 5?. is the z-component of the spin op- 
erator at site i and imaginary time r. In Fig. the 
distribution function of is shown. The distribution 
function has two peaks; one sharp peak at —b and the 
other broader peak at b/2. This is what we expect for 
the VBS-H phase. For the frozen plaquette phase, one 
expects three peaks at a, 0, and —a. 

The VBS-H phase also has magnetic order. In Fig. 21 
we show the finite size scaling of the equal time spin-spin 
structure factor at q = Qo, 

S{q)=N{SlM^ (4) 

where ^ = (V^) ^fr exp(^qr,). 

Therefore, our Monte Carlo data are consistent with 



FIG. 7: (color online). The correlation function Cb(ri — r^) 
between the bond indicated by a large text (red online) and 
the other bonds for L — 12, Jz/J± = 10.5, D = 0, and 
T = 0.05Jx/. Resonating triangles (trimers) are denoted by 
thick lines. 

the resonating plaquette phase (VBS-H) . It is worth not- 
ing that the VBS-H phase is quite similar to the VBS 
phase discovered in the hard-core boson model on the 
kagome lattice at fillings of 1/3 and 2/3 (for details, see 

B. VBS-T(B) phase 

Another interesting result in this work is the discovery 
of another phase with resonating plaquettes at small val- 
ues of D/ Jx and large values of Jz/J± as shown in Fig. [T] 
This phase was not predicted before. There are actually 
two different competing phases that are very close in en- 
ergy. Triangles resonate in one of those phases (VBS-T 
phase) and bow-ties resonate in the other one (VBS-B 
phase). We are not able to determine reliably which 
phase is the true ground state. 

We find that there is a finite temperature transition 
from the high temperature paramagnetic phase to the 
resonating phase for J^/ J± > 10 and D close to zero. We 
have not attempted to obtain the precise location of this 
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FIG. 8: (color online). The correlation function Cb(ri — rs) 
between the bond indicated by a large text (red online) and 
the other bonds for L — 12, Jz/J± ~ 10.5, D = 0, and 
T — 0.002Jx. Resonating bow ties are denoted by thick lines. 



transition. The transition temperature Tq ~ J±/10 for 
JJJ± = 10.5 and D = 0. Typically, we find the VBS-T 
phase just below the transition and either the VBS-T or 
VBS-B phase at much lower temperatures depending on 
the configuration we start our Monte Carlo simulations 
with. 

The bond-bond correlators have well pronounced peaks 
at Qi = (tt, 0), Qa = (0, tt), and Q3 = (tt, tt) in the VBS- 
T phase and at those or other points in the VBS-B phase. 
In Fig. [HI the finite size scaling of the equal time bond- 
bond structure factor S'b(Q3) defined in Eq. [3]is shown 
for the VBS-T phase. The structure factor divided by 
the number of sites scales to a finite value in the ther- 
modynamic limit indicating long range valence bond or- 
der. Note that the VBS order parameter is very large^^ 
in sharp contrast to the VBS-H bond order parameter 
and to some other models, in which the VBS order was 
confirmed in quantum Monte Carlo simulations.^^ To in- 
vestigate the nature of the VBS-T(B) phase, we study 
the real space correlation function defined in Eq. [2l As 
shown in Figs. [7] and [51 the VBS-T phase exhibits a net- 
work of resonating triangles and the VBS-B phase shows 
a network of resonating bow-ties. 

It is worth noting that the presence of the insulat- 
ing phase at _D = is in sharp contrast to the spin- 1/2 
XXZ model with ferromagnetic XY and antiferromag- 
netic Ising exchange interactions on the kagome lattice, 
where the uniform XY ferromagnet persists at any finite 
value of Jx/J^Jiii^ 

The appearance of the VBS-T or VBS-B states can be 
qualitatively explained as follows. Let us first consider 
the D — case. The classical model {J± = 0) has the 
ground state constraint such that each triangle on the 
kagome lattice should have one of the spin arrangements 
++-, + — , or +0- (and permutations), where the eigenval- 
ues of the spin-1 operator are denoted by +, -, and 0. 
There are many different ways to arrange classical spins 




FIG. 9: (color online). A sample configuration in which 
the central triangle can resonate without leaving the classical 
ground state manifold. 



on the lattice to fulfill the above constraint, leading to a 
highly degenerate classical ground state manifold. The 
third group of states (+0- and permutations) is quite 
special in the quantum case — a kinetic term (or an XY 
term) acting on such a state leaves the spins on a triangle 
(a flippable triangle) in the classical ground state mani- 
fold; for example, 5'^S';7|0+~) = I+0-). In general, such 
a move may violate the classical ground state constraint 
on the neighboring triangles. However, it does not vi- 
olate the constraint in the special case when the three 
neighboring triangles of a flippable triangle have + and 
- spins as shown in Fig. [9l Thus the flippable triangles 
can resonate without leaving the classical ground state 
manifold. This consideration naturally leads to the fol- 
lowing wave function that has a kinetic energy gain of J± 
(over the classical ground states) in first order (degener- 
ate) perturbation theory 

ItP) = |0+-) -I- |0-+) + I+0-) + |-0+) -I- I+-0) -I- |0-+). 

On the other hand, two flippable triangles shown in 
Fig. [9| can resonate simultaneously giving rise to a res- 
onating bow-tie. It is expected that turning on an in- 
finitesimal J± will favor the spin configurations where 
either the number of independently resonating triangles 
or the total number of resonating triangles is maximized. 
In the latter case, we expect that resonating bow-ties 
are close-packed as shown in Fig. [5) This arrangement 
of bow-ties violates the classical ground state constraint 
and one needs to project out the components of the wave 
function that violate the constraints^ It is not possible 
to determine which state has lower energy based on the 
above analysis. 

Therefore, quantum fiuctuations lift the macroscopic 
degeneracy of the classical ground states and gives rise to 
a VBS phase via quantum order by disorder mechanism. 
The VBS phase should survive at small but finite values 
oi D/Jz because of a finite excitation gap. The state with 
maximal number of independently flippable triangles is 
exactly the VBS-T state that is found in quantum Monte 
Carlo simulations, i.e. the configuration shown in Fig. [T] 
The spin configuration, where all the flippable triangles 
resonate, is the VBS-B phase and is shown in Fig. [5] The 
energies of those two states are very close. 
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FIG. 10: (color online). Finite size scaling of the superfluid 
density ps at Jz/J±_ = 6 for jS/L = 1/Jx- 



4 

3.5 
3 

2.5 
2 

1.5 
1 

0.5 





-0.5 







0.5 

l/v 



1.5 



FIG. 11: (color online). Data collapse of the superfluid den- 
sity Pa at 3^1 3i_ = 6 for /3/L — 1/J_i. 



As we have noted above, the bond order is very strong 
in the VBS-T(B) phase. This may be explained by the 
fact that the VBS-T phase can be chosen via resonat- 
ing triangles in first order perturbation theory whereas, 
in most of other cases (including VBS-H), plaquette res- 
onance is obtained in higher order perturbation theory. 
As a result, the VBS-T(B) phase may be more robust 
than other cases. 



C. Superfiuid-paramagnet phase transition 

In this subsection, we discuss the phase transition be- 
tween the superfluid (or XY ferromagnet in spin lan- 
guage) and the featureless quantum paramagnet. Our 
numerics shows that this transition is continuous and 
most likely belongs to the 3d XY universality class. We 
analyze the finite size scaling of our data as follows. We 
measure the superfiuid density by measuring the wind- 
ing number fluctuations.^^ In the vicinity of a continuous 
transition, the superfluid density scales as 



(5) 



where Fp^ is the scaling function, L is the linear system 
size, z is the dynamical critical exponent, v is the corre- 
lation length exponent, — K = [D / Jz)c — D / Jz \s the 
distance to the critical point, and (3 is the inverse tem- 
perature. To cross the phase boundary, we change D and 



FIG. 12: (color online). The superfluid density ps as a func- 
tion of Jz/J_L for L = 12, D = 0, and different temperatures. 
Lines are guides to the eye. 



keep Jz fixed. The data scale very well with the dynami- 
cal critical exponent z = 1. In Fig. 1101 the superfluid den- 
sity Pa times the system size L is shown as a function of 
the coupling constant. As follows from the above scaling 
form, the curves for different system sizes should cross at 
the transition point if the ratio is fixed. Such a dis- 
tinct crossing point is seen at [D/ Jz)c — 1.2124. It also 
follows from Eq. [5] that the curves for different system 
sizes should collapse onto a universal curve for appropri- 
ate values of v and {D/ Jz)c when psL is plotted as a func- 
tion of [{D/Jz)c - D/Jz\L'^'''. In Fig.ini we show such 
a data collapse for v = 0.67(2) and {D/Jz)c = 1.2124(2). 
The error bars are estimated from the stability of the data 
collapse with respect to varying the fitting parameters. 
Thus one may conclude that the superfiuid-paramagnet 
quantum phase transition is continuous. One may also 
infer that the transition is in the 3d XY universality class 
from V = 0.67(2). However, we have not measured the 
other critical exponents that could be used to unambigu- 
ously confirm this prediction. 



D. Superfluid- VBS phase transitions 

The quantum phase transitions out of the superfluid 
phase to VBS phases are strongly first order. As an 
example, let us consider the transition from the super- 
fiuid phase to the VBS-T phase. In Fig. [T^l we show 
the superfluid density ps as a function of Jz/ Ji_ at dif- 
ferent temperatures. The superfluid density jumps at 
large values of Jz/ Jjl indicating a transition to an insu- 
lating phase (VBS-T phase). We confirm that the transi- 
tion is strongly first order by observing hysteresis effects 
upon increasing or decreasing Jz/ J_l across the transition 
(not shown) and by a double peaked structure in the dis- 
tribution of the XY energy (kinetic energy of bosons), 
J±{SfSj + SfSj), close to the transition (not shown). 

The transition from the superfluid phase to the VBS-H 
phase is also first order. In the top part of the VBS-H 
lobe (see Fig. [1]), the flrst order nature becomes some- 
what weaker as the value of J±/Jz is decreased. We do 
not observe a direct transition from the VBS-H phase to 
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the quantum paramagnetic phase, which was predicted 
in Ref. [HI, in the parameter region we can access, i.e. for 
J±/Jz <i 0.08. There is a narrow region of the super- 
fluid phase between those two phases, as shown in Fig.[T] 
However, we cannot rule out the possibility that there 
might be a direct transition at much smaller values of 

III. DISCUSSION 

We now discuss possible origins of the discrepancy be- 
tween our quantum Monte Carlo results and those of 
previous analytic approaches with various approximation 
schemesj^i^i^ In particular, it is found that the VBS-H 
phase, not the FP, is the stable ground state for mod- 
erate strength of D/J^ and small J±/Jz- This suggests 
that the string-net mean field theory analysis in a previ- 
ous work^ may not be sufflcient for the identification of 
the true ground state. 

The string-net picture starts from an alternative rep- 
resentation of the spins on the kagome lattice, where the 
spins are considered to be in the middle of the links con- 
necting the sites of the honeycomb lattice (that can be 
obtained by connecting the centers of triangles on the 
kagome lattice).- This honeycomb lattice is a bipartite 
lattice and consists of A and B sublattices. The occu- 
pation of a given link / = {ij) by a string is defined as 
follows; the link contains an oriented string pointing from 
i e A to j e i? if 5| = +1 and from j to i if 5| = -1. 
The link is empty if 5| = 0. Then the spin configura- 
tions on the original kagome lattice and those of oriented 
closed strings are in exact correspondence.^ 

The string-net mean-field theory uses the following 
ansatz for the variational ground state wavefunction.^ 

where X represents an oriented string configuration and 
{zij} a large number of variational parameters, and riij 



the occupation number of the oriented link ij. The mean- 
field analysis begins with a (translationally-invariant) 
string liquid state (or a string condensed state) where 
can be set to a constant a. Considering the spectrum of 
collective modes in this state and the identification of the 
wavevector where the collective modes become soft, one 
can in principle study the instability to a translational- 
symnietry broken state. More explicitly, this can be 
achieved by writing = ae^*^^*'^'^ and studying the 
fluctuations of Eij and Aij. It was found that the soft 
mode is described hy E ^ ^E+q + ^*E-q and A = 0, 
where Eq is the Fourier mode at Q = (47r/3,0) and 
equivalent wavevectors.-- Here $ is a complex number. 

The energy (or the Ginzburg-Landau theory) of the 
system as a function of <& can be obtained as^ 

i7($) = + + C[$^ + ($*)^] + • • ■ , (7) 

where A, B, C are real constants and (f> can be regarded 
as an order parameter. The choice of the ground state 
is sensitive to the phase of $ and hence the sign of C 
basically determines the true ground state. It was found 
that if C is positive (negative), then the FP (VBS-H) is 
favoredi^ This sign, however, is very difficult to deter- 
mine in analytic approaches. In Ref. 0, an unrestricted 
variational wavefunction calculation was also done on a 
small system size 3x3, leading to the conclusion that the 
ground state may be the FP phasci^ We think, however, 
that this system size is perhaps too small for definitive 
conclusion. Indeed our quantum Monte Carlo results are 
clearly consistent with the VBS-H phase, and hence the 
negative sign of C in the Ginzburg-Landau theory. 
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